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Abstract Compacted unbound granular materials are ex- 
tensively used as sub-layer in pavement design. Most pave- 
ment design guides assume that they are responsible for 
the degradation and deformation of the roads and rail- 
ways that they support. Biaxial tests are usually employed 
to investigate the elasto-plastic response of these materi- 
als to cyclic loading. A particularly interesting question 
is whether a limit load exists, below which the excita- 
tions shake down, in the sense that the material does not 
accumulate further deformations. We have carried out a 
detailed study of the elasto-plastic behavior of a simple 
model of unbound granular matter submitted to cyclic 
loading. The dissipated energy through out the simula- 
tion has been used for the characterization of the different 
regimes of responses. 
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Introduction. 

Traditional pavement design methods are still almost com- 
pletely empirical. Long term experience with the perfor- 
mance of on-service roads is supplemented with the re- 
sults obtained from especially constructed test pavements. 
Changes in material properties, the loading or the envi- 
ronmental conditions can be hardly investigated. The in- 
formation obtained from the experiments is therefore lim- 
ited. The disadvantages of traditional design have become 
more obvious during the last decades as a consequence of 
the growth of transportation needs and the urge for us- 
ing recycled materials. Mechanistic and analytical design 
procedures have consequently been developed based on 
the analysis of the response of the structure under specific 
loads. Understanding the behavior of the components of 
the structure is of course a pre-requisite of this approach. 
In this respect, the deformation behavior of unbound gran- 
ular materials (basic component of roads and pavements) 
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has been one of the main topics of pavement engineering 
for many years [1-3], since they are principal responsible 
of the rutting and cracking of the pavement [4] . 

Grain scale investigation of cyclic loading of granu- 
lar materials is possible using discrete element methods, 
that solve the dynamical evolution of the system accord- 
ing to the interaction between the particles. The quasi- 
static Coulomb's friction on lasting contacts has been ex- 
tensively studied by means of algorithms of Molecular Dy- 
namics (MD) [5]. In this method, a visco-elastic interac- 
tion is introduced in each contact. The evolution of the 
grains is solved explicitly by the numerical solution of the 
equations of motion. The grains are usually represented 
by spheres (in 3D) or disks (in 2D). One expects that 
this approximation will reproduce the behavior of smooth- 
grained unbound granular material (UGM) , such as grav- 
els. Materials like crushed aggregates, composed by grains 
with more irregular shapes, have also been modeled [6-9], 
leading to more complicated and therefore less efficient 
algorithms. 

It is usually assumed that the global response of an 
UGM to cyclic loading can be decomposed into resilient 
deformation, and permanent deformation which eventu- 
ally leads to an incremental collapse of the pavement. Note 
that even in the case of very small permanent deforma- 
tions appearing in each cycle, its systematic accumulation 
could lead to an eventual failure of the structure due to 
excessive rutting. Whether a given material will experi- 
ence progressive accumulation of permanent deformation, 
or whether this process will stop is therefore crucial for 
the performance prediction. There arc experimental evi- 
dences pointing out that permanent deformations quickly 
shake down, becoming the response basically resilient after 
a certain adaptation period [10,11]. A micro-mechanical 
observation of the permanent and resilient deformation 
will help to get insight into the mechanisms involved in 
the shakedown. 

Most of the theoretical approaches to UGM deforma- 
tion try to identify the internal variables of postulated 
constitutive equations of the material, based on macro- 
mechanical observations of the response of soil samples 
in triaxial or biaxial tests. Shakedown theory, however, is 
concerned with the evolution of the plastic deformation in 
the material. It predicts that a structure is liable to show 
progressive accumulation of plastic strains under repeated 
loading if the magnitude of the applied loads exceeds a cer- 
tain limiting value, the so-called shakedown limit or limit 
load. The structure is then said to exhibit an incremental 
collapse. On the other hand, if the loads remain below this 
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Fig. 1. Visco-elastic stress-strain cycle. The resilient and per- 
manent deformation is clearly identified. At the end of the 
cycle, the system did not return to its original position, but 
accumulated a certain permanent deformation e. The energy 
dissipated in each cycle is the area enclosed by the loading and 
unloading curves. 



limit, the growth of plastic deformations will eventually 
level off and the structure is said to have attained a state 
of shakedown by means of adaptation to the applied loads. 
Under these premises, four categories of material response 
under repeated loading can be distinguished [12]: 

— Elastic range for low enough loading levels, in which 
no permanent strain accumulation occurs. 

— Elastic shakedown. The applied stress is slightly below 
the plastic shakedown limit. The material response is 
plastic for a finite number of cycles, although the ulti- 
mate response is elastic. 

— Plastic shakedown. The applied stress is low enough 
to avoid a rapid incremental collapse. The material 
achieves a long-term steady state response with no 
accumulation of plastic strain, but hysteresis in the 
stress-strain plot. 

— Incremental collapse. The repeated stress is relatively 
large, so that plastic strain accumulates rapidly with 
failure occurring in the relatively short term. 

The aim of this paper is to prove the existence of 
a shakedown regime in a simple model of UGM. This 
is done characterizing the material response to different 
loadings by means of the calculation of the dissipated en- 
ergy throughout the experiment. A further description of 
the shakedown in the framework of the general elasto- 
plastic theory will be presented in the next section. The 
model used in our simulations is carefully described in 
Section 3, where the inter-particle interactions are also 
detailed. The results of the simulation of the system sub- 
mitted to different loadings ranges are shown in Sections 
4 and Section 5 . Finally, in Section 6 the conclusions of 
this work arc summarized. 



Shakedown theory. 

Shakedown analysis is basically an extension of limit anal- 
ysis to the case of periodic loading [13]. It was introduced 
in the context of structural design to take into account 
that, in the case of repeated loads, an accumulation of 
small plastic deformations may occur in every cycle, that 



eventually might lead to a collapse of the structure. It 
was necessary, therefore, to develop methods able to es- 
timate the limit load, beyond which failure occurs by ei- 
ther incremental collapse or fatigue. Below this load the 
structure will shake down, in the sense that plastic de- 
formation ceases at some point or, more precisely, that 
the energy dissipated is bounded in time. Melan's pioneer 
works [14] were completed by Koiter [15], who first formu- 
lated a general kinematical shakedown theorem. Melan's 
and Koiter's theorems offer an lower and upper bond for 
the limit load, respectively. Their ideas are applied in [3] 
to a simple pavement model subjected to repeated moving 
loads, transforming the shakedown theorems into a linear 
programming problem, which is then numerically solved. 

Shakedown theory should be understood within the 
general framework of the classical theory of elastoplas- 
ticity, where the existence of an elastic regime is postu- 
lated [16]. Elastoplastic constitutive laws consist of linear 
strain-stress relations that result in a overall incrementally 
non- linear relation. The basic assumption of these theories 
is that the stress state at any point of the material can 
be separated into two clearly different elastic and plastic 
parts: 



i(x,y) = o-fJx,y) +a[Ax,y), 



(1) 



where afj represents the residual stress due to the accu- 
mulated permanent deformation (see Figure 1) and afj is 

the elastic stress, related to the resilient strain through 
the generalized Young's moduli E^i : 

= (2) 
This relation is usually known as the generalized Hooke's 
law. In the space of stresses, there is a certain volume 
around the origin in which the system experiences an elas- 
tic deformation. The surface enclosing this area is a very 
special surface, called the yield surface. It is calculated in 
terms of the yield function f((Tij, efj): 

/K-,e£)=0, (3) 

The yield function / is characteristic of the material. If it 
does not depend on the accumulated deformation, , the 
material is ideally plastic, whereas in any other case, it is 
said to suffer hardening. As hardening occurs, the size of 
the yield surface consequently changes, but in the case of 
isotropic hardening, neither the shape nor the orientation 
of the yield surface is altered. 

When the stress state reaches the yield criterion (3) 
the material undergoes plastic deformations, being the in- 
crements of plastic strain given by the normality rule: 

< = <-> 

In this expression, the existence of the plastic potential 
function Q, to which the incremental strain vectors are 
orthogonal, has been assumed. A<P determines the magni- 
tude of the plastic deformation. In the simplest case, the 
plastic potential Q and the yield function are the same, 
i.e. Q = / (associative flow rule of plasticity). It will be 
assumed in the following that, during the plastic deforma- 
tions, the yield surface translates in the stress space. This 
can be expressed: 

-aij,efj) = 0, (5) 
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where a, 3 represents the total translation of the yield sur- 
face. Note that (4) and (5) are special cases of general 
anisotropic hardening, involving change in size, shape and 
orientation of the yield surface. Let us now consider the 
experiment shown in Figure 1, in which the system reacts 
plastically to the imposed load. The response of the sys- 
tem is then discretized in r linear segments in which the 
system behaves as described in equations (4) and (5) . The 
total strain can still be decomposed in a plastic and an 
elastic component, as in (1), and the plastic strain rate 
can be determined, 



fc=i 



13 ^ d*v' 



(6) 



in terms of the r yield surfaces fk and the material rates 
tf rfc , which respectively fulfill the following conditions: 

i> k >0, & k f k = forfc= l,...,r. (7) 

The yield condition has been then separated in r parts, 
each of them given by: 



fk(vn-a!ij) < 0. 



(8) 



Physically, it has been assumed that there exists in the 
domain of interest a series of yield surfaces, each defining a 
region in the total stress space of the system. Note that, as 
a consequence of this construction, the system will behave 
elastically if it is subjected to a low enough stress. Note 
also that surface /o is contained in surface fi, which is 
contained in /2, and so on... . In the experiments, the 
surfaces may deform and translate in the stress space. As 
soon as surface fk touches fk+i, they will move together 
until they touch the next surface. The shift of fk, ot k j, will 
be proportional to the plastic displacement in the set of 
yield surfaces, and the same reasoning can be applied to 
its rate, 



i=i 



dfl_ 
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(9) 



being c kl the strain-hardening coefficients. For this model, 
under the proper conditions [13], Melan's static shake- 
down theorem assures that the total energy dissipated in 
any allowable stress path is bounded [12]. Since direct ex- 
perimental measure of the dissipated energy in the cycle 
is hardly possible, computer simulation seems to be the 
best tool for investigating the shakedown. In the particular 
case of unbound granular materials, Molecular Dynamics 
is the more appropriate method, for it solves the dynam- 
ics of the system after well established collision rules, and 
also because a detailed study of the micro-mechanics of 
the system and its relation to the macroscopic behavior of 
the material is possible. 

Werkmeister et al. [17] have presented some exper- 
imental results on the shakedown of unbound granular 
material subjected to triaxial tests. They conclude that 
shakedown theory can be applied to study the formation 
of instabilities and degradation of UGM, only if some 
aspects of it are modified. The existence of three pos- 
sible responses of the system is consequently proposed, 
depending on the load. For low loads, the system will 
definitely shakedown (behaving almost elastically), after 



a post-compaction transient in which plastic deformation 
is accumulated. On the contrary, for very high loads the 
material deforms boundless at an almost constant rate. 
The case of intermediate loads is somehow ambiguously 
treated by Werkmeister. Experimentally, however, it is 
clearly found that the plastic strain rate decreases to a 
low, nearly constant level for intermediate loads. A small 
residual linear increment of plastic strain is accumulated 
after each cycle. 



3 

Model. 

A polydisperse system of inelastic disks has been used to 
model the UGM. The chosen type of grain is an appropri- 
ate idealization of natural sands like gravel, in which the 
aggregate particles are smooth and round. The boundary 
conditions of the biaxial test will be reproduced. The sys- 
tem is bounded by four mobile walls that induce a certain 
stress state in the system. The grains interact with each 
other through a viscoelastic force that will be explained 
in detail next. The dynamics of the system is solved by 
means of a MD algorithm. 



3.1 

Contact forces. 

In order to calculate the forces, it is assumed that all the 
disks have a characteristic thickness L. The force between 
two disks is then written as F = fL and the mass of 
the disks as M = mL. In any contact between two real 
grains, there is deformation in the impact region. In the 
simulation presented here the disks are supposed to be 
rigid, but they can overlap so that the force is calculated 
by means of this virtual overlap. 

The contact force can be divided in two components, 
f c = f e + f v , where f e and /" are the elastic and viscous 
contribution. The elastic part of the contact force also 
splits in two components, / e = f^rf + fjrt c . Where n c = 
(rj — rj)/\ri — Tj\ is the unitary vector pointing from the 
center of mass of particle j to particle j and t c — u z x h c 
(u z is the unitary vector in the direction perpendicular to 
the plane of motion of the disks). 



3.1.1 

Normal elastic forces. 

Given two overlapping disks i and j whose diameters are 
di and dj , respectively, the normal elastic force acting on 

Where k n is the normal stiff- 



1 di+dj 



them is f^ = 

ness and A is the overlapping area. The overlapping area 
of these two particles can be easily calculated in terms 
of their diameters and the distance that separates their 
centers, ya; 



A 



-2y x + d 2 j 




+ arcsin (xo/dj) 



arcsin (xo/di) 



+ 



(IP) 
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where the intersection point of both circles, x , has been 
introduced, 



2y 



(ii) 



Note that the use of this model for the calculation of the 
elastic force is similar to the use of the overlapping length 
[18], except for a renormalization factor in the clastic con- 
stant of the spring. 



3.1.2 

Friction forces. 

In order to model the quasi-static friction force, the elas- 
tic tangential force is calculated using an extension of the 
method proposed by Cundall-Strack [6]: An elastic force 



ft 



-k t Ax t proportional to the elastic displacement is 



included at each contact. k t is the tangential stiffness of 
the contacts. The elastic displacement Axt is calculated 
as the time integral of the tangential velocity of the con- 
tact during the time that the elastic condition |/ t | < \ij n 
is satisfied. Here p is the friction coefficient. The slid- 
ing condition is imposed keeping this force constant when 
| ft | = fifn . The straightforward calculation of this elastic 
displacement is given by the time integral starting at the 
beginning of the contact 



Ax e t 



(12) 



where O is the Heaviside step function and v% denotes the 
tangential component of the relative velocity v c at the 
contact. v c depends on the linear velocity Vi and angular 
velocity uji of the particles in contact according to: 



V = V; — Vj — U), 



di „ c 
x-n 



+ u}, x 



do 



(13) 



3.1.3 

Damping forces. 

A rapid relaxation during the preparation of the sample 
and a reduction of the acoustic waves produced during 
the loading is obtained if damping forces are included. 
These forces are calculated as /" = —vvnv c . Here m = 
mivrijl (nii + nij) is the relative mass of the disks in con- 
tact and v is the coefficient of viscosity. These forces in- 
troduce time dependence effects during the cyclic loading. 
Nevertheless, these effects can be arbitrarily reduced by 
increasing the time of loading, so that the quasi-static ap- 
proximation can be assumed. 



3.2 

Sample preparation. 

The disks are initially randomly distributed into a rect- 
angular box which is big enough for them not to overlap. 
Their radii distribution is Gaussian with mean value A 
and standard deviation 0.23A. The interaction of the disks 
with the walls of the box is implemented applying a nor- 
mal visco-elastic force f™ = —k n S — m^„v£ at each disk 
lodged in any of the walls. Here 5 is the penetration length 




Fig. 2. Hambly's principle for biaxial test. The degrees of 
freedom of the walls allow to impose any pair of stresses <7i, 
<J2 to the system. 



of the disk into the wall. A gravitational field g — gr is also 
temporarily switched on, where r is the vector connecting 
the center of mass of the assembly with the center of the 
disk. The activation of this gravity field produces homoge- 
neous, isotropic distribution of disks. After a certain time 
to, (which is defined below) a modulation in gravity is im- 
posed such as g = g + l/2(gj — <7o)(l + cos (WOnt / to)) 
until the time 2t in order to reduce the porosity. For 
g/ = lOOgo (go will be defined next), samples arc obtained 
with packing fraction 0.841 ± 0.001. 

The external load is imposed applying a force a x W and 
<7 y H on the horizontal and vertical walls, respectively. W 
and H are correspondingly the width and the height of 
the sample. When the velocity of the disks vanishes, the 
gravity is switched off. A fifth-order predictor-corrector 
algorithm is used to solve the equations of motion. The 
time scales used both for the compression and the cyclic 
loading are such that we can assume a quasi-static regime. 
This allows us to consider a static friction coefficient only. 

It can be shown by dimensional analysis that the strain 
response depends only on a minimum set of dimension- 
less parameters: 1) the ratio to As, between the period of 
cyclic loading t and the characteristic period of oscilla- 
tion t s = \Jknl ' p\ 2 (where k n is the normal stiffness of the 
contacts and p the density of the grains). 2) The ratio t r /t s 
between the relaxation time t r = 1/u and the oscillation 
time. 3) the ratio k t /k n between the stiffnesses. 4) the adi- 
mensional stress state a/k n and 5) the friction coefficient 
\x. In our simulation the following values have been kept 
constant: to = 1000t s , t r = lOOt^, k t — 0.33fc„, k n = 2 
MPa, \i = 0.1, go = 6.25 x 10~ 8 fc„, the initial pressure 
Po = 6 x 10- 4 fc„, and the MD time step t s = 2.5 x 10~ 6 s. 
The imposed load, Aa, has been the variable parameter 
used to investigate the response of the system. 



3.3 

Biaxial test. 

The boundary conditions in the simulations are those of a 
biaxial test according to Hambly's principle (see Figure 2). 
In this experiment, the probe is compressed by four mobile 
walls which exert on the system a fixed stress o\ = a xx in 
the x direction and oi = a yy in the y direction. The stress 
values can be changed and the corresponding changes in 
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the principal components of the strain tensor ei, e 2 can be 
measured. In our tests, the system is first homogeneously 
compressed with <j\ = a 2 . After the system has reached 
an equilibrium state under the pressure Po = <Jl +' J2 = <j\ , 
the vertical stress is quasi-statically changed: 

Aa / (2wC 

— — 1 — cos 

2 V V t 

where t is the simulational time and to is the period of the 
cyclic loading. It will be useful to define the deviatoric 
stress of this configuration, SV,-: 



(J2{t) = Po 



1 + 



(14) 



Si(t) 
S 2 (t) 



<Tl(t)-<T2(t) 









PoAa 

4 



(l-cos(^)) 




PpAcr 



(l-cos(^)) 



(15) 



In the following experiments, the changes in the load- 
ing will be characterized by the parameter Aa, which is 
directly related to the maximum value of the deviatoric 
strain. The response of the system will be characterized 
by the adimensional quantity 7, which is defined in terms 
of the deviatoric permanent strain. This, analogously to 
the deviatoric stress, is the difference of the strains in the 
principal directions. The permanent strains in the princi- 
pal directions are: 

L x (t) 



<*(*) - 



LI ' 

Ly{t) 



(16) 
(17) 



Where L x / y (t) are the dimensions of the system at the 
time t, whereas L^/y are the dimensions at the beginning 
of the cyclic loading. Then, 7 is defined as: 

7 = e 2 -ei. (18) 
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Fig. 3. Dissipated energy per cycle A in the system for differ- 
ent loadings Aa. After a transient in which energy is quickly 
dissipated, the system reaches a state in which A is approxi- 
mately constant. This value is higher for higher Aa. There is 
however a critical value of Aa below which there is only viscous 
dissipation in the system. 



Description of the plastic shakedown. 

According to shakedown theory, there are two possible 
fully clastic ranges in the deformation of the material. On 
one hand, a pure clastic response will be found for low 
enough loadings. On the other hand, clastic shakedown 
will come up if the load is slightly increased: The system 
will show some plastic deformation, that will eventually 
level off leading finally to an elastic response. 

The response of the system in the limit Aa — > will be 
discussed at this point. Within the framework of elasto- 
plasticity, the material should then behave elastically, be- 
ing the recoverable strain after a cycle a mathematical 
function of the stress. In the simplest case this relation 
should be linear, Eq.(2). As pointed out in Figure 1, the 
dissipated energy is the area of the cycle (C) in the stress- 
strain plane. Note that the total work applied to the sys- 
tem is the integral of the stress-strain curve during the 
loading (L). The ratio of these two magnitudes will help 
us to characterize the response of the system: 



A = 



§c g(e ) & 

S L a ^) de 



(19) 



Figure 3 shows the evolution of A on time as Aa is 
changed. The system dissipates the most energy in the 
first cycles, reaching a state in which the energy dissi- 
pated in each cycle remains constant. This reflects the 
similarity of the stress-strain cycles and is associated to 
a certain periodicity of the sliding contacts [7]. The final 
value of A clearly depends on Aa. Note also that it takes 
longer for the system to reach this final stage if the load- 
ing is increased. In fact, if the loading is very low, the 
relaxation is apparently immediate. In that situation, the 
limiting value reached by A is independent of the load- 
ing imposed. For Aa = 10~ 3 and Aa = 10~ 4 , the energy 
dissipated in the system converges to the curves corre- 
sponding to Aa < 10~ 5 . This decay is associated with the 
disappearance of sliding contacts in the sample. These re- 
sults confirm the existence of a range of possible loadings 
for which the response of the UGM would be a plastic 
shakedown. The system has accumulated systematic plas- 
tic deformation at the beginning of the loading, but this 
process stops when the contacts cease sliding. In the limit 
of Aa — > one would expect that the dissipated energy 
per cycle also tends to zero. This is in fact what is found 
in the simulation, but Figure 3 also shows that the dissi- 
pated energy tends to zero keeping the ratio A constant. 
This is associated to the fact that no sliding contacts are 
found in the system for that range of excitations. The dis- 
sipation is therefore completely due to the viscosity. The 
energy dissipated through the viscosity in this model is 
directly related to the mean velocity of the particles. In 
the ideal quasi-static limit, this velocity would be zero. In 
our simulations, it can be reduced making to bigger. For 
a finite to, however, no purely elastic regime can be found 
in this model. 

Let us next investigate the deformation of the system 
in the viscous regime. Part (a) of Figure 4 shows the evolu- 
tion of the permanent strain in four different experiments 
all of them corresponding to very low load. In the lowest 
load case, Aa = 10~ 8 , the permanent strain does not seem 
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Fig. 4. (a) Evolution of the permanent plastic strain in sim- 
ulations with different deviatoric strains. For comparison, the 
values of the strain are scaled with the initial value 70. (b) 
Stress-strain plot for the case Aa = 1CP 8 . The dotted line is 
plotted only as reference, showing a linear behavior. The de- 
tails of the simulations are given in the text. 

to have a defined trend, but on the contrary, it apparently 
wanders around the initial value, reflecting the existence 
of dissipation (Part (b) of Figure 4). The deviations from 
the original value grow as Aa increases. For the highest 
load shown, a steady expansion of the system follows a 
small initial compression. In the intermediate cases, the 
strain remains approximately constant until N = 500 in 
the case Aa = 10~ 7 , whereas for Aa = 10~ 6 , the change 
appears even earlier. In part (b) of the figure, the stress- 
strain cycle is plotted for the lowest loading and compared 
with a linear function. The existence of dissipation can be 
observed, although the ratio of dissipated energy and the 
total elastic work is A = 0.017. 



Investigation of the plastic response. 

The plastic deformation induced in the granular mate- 
rial due to the cyclic loading has also been investigated. 
Given a system configuration and a certain confining pres- 
sure Pq the response of the system to different deviatoric 
stresses is shown plotting the strain rate Aj/AN versus 
the permanent strain 7, Figure 5. Part (a) of the figure 



helps to identify easily the existence of incremental col- 
lapse in the sample. For high loads, the strain rate does 
not diminish after each cycle, but remains approximately 
constant or even grows, reflecting the continuous accumu- 
lation of plastic deformation that the system is bearing. 
For low loading, however, the strain rate quickly decays 
to a certain almost constant value, implying a linear in- 
crease of the strain with the number of cycles. In both 
cases, a final constant value of the strain rate can be mea- 
sured, A^p/AN. Note that the existence of a shakedown 
of the excitations cannot be consequently inferred from 
this graph. On the contrary, this ratcheting regime is as- 
sociated with a periodic behavior of the sliding contacts 
and has been already reported both in a MD simulation of 
a polygonal packing [7] and experimentally [17]. The final 
permanent deformation is higher for higher applied loads. 
This is related to the dependence of the strain rate on Aa 
(Figure 5 (b)). The strain rate is higher as the loading is 
increased. The incremental collapse is clearly associated 
to an unbound growth of the strain rate. 

These results are compatible with experimental mea- 
surements on a triaxial test by Werkmeister et al. [17], 
even though no grain crushing is included in our model. 
This is probably the reason why the incremental collapse 
seems to be weaker in the simulation as it is in the exper- 
iment. 

The evolution of the strain rate for intermediate load- 
ings is worth further discussion. Figure 6 shows the re- 
sponse of the system for different excitations below the 
incremental collapse. Let us remark first the existence 
of two different regimes in the response of the material. 
At the beginning of the loading process, the system re- 
acts accumulating deformation at a (relatively) high rate. 
This regime will be called post- compaction. After the post- 
compaction, the material has adapted to the new situation 
and the strain rate is very low and almost constant. This 
final stage is the ratcheting regime, in which the mate- 
rial accumulates deformation linearly. Within the loading 
cycles the system behaves quasi-periodically. The strain 
growth is linked to a small compression of the system. As 
Aa is increased, the duration of the post-compaction (in 
terms of the number of cycles) is longer. This is related to 
the existence of transient states, in which the strain accu- 
mulates at a constant rate for some cycles, before reaching 
its final level. This intermediate state is clearly observed 
in cases Aa = 0.31, 0.32 and 0.34. The strain rate decays 
to a fixed value for some time, but then the material reacts 
and an abrupt change on the strain rate leads to the final 
relaxation. This phenomenon is less clear as Aa increases 
and the strain is accumulated at an almost constant rate 
during the post-compaction. In this situation, if Aa is big 
enough, the system can even break, as it is noticeable in 
the case Aa = 0.36 in the figure. During the experiment, 
a sudden compression of the system occurs just after the 
post compaction, and previous to the complete adaptation 
of the material. 

Extending the results of Figure 3, a similar study on 
the dissipated energy in shown in Figure 7 for higher load- 
ings. According to the shakedown theorems, none of the 
material to which these curves correspond would shake 
down, for the total dissipated energy is not bounded. The 
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Fig. 5. (a) Strain rate versus the accumulated permanent 
strain for different Aa. (b) Dependence of the final value of 
the strain rate (Ajf / AN) on the loading Aa. All the simula- 
tions correspond to the same initial configuration of 400 disks. 



Fig. 7. (a) Evolution of the scaled dissipated energy A in time 
for different values of Aa. The plotted values correspond to 
(from bottom to top) Aa = 0.01, Aa = 0.10, Aa = 0.20, 
Aa — 0.35, and Aa — 0.45. (b) Final dissipated energy ratios. 
Simulation details are the same as in Figure 5. 
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Fig. 6. Decrease of the strain rate through the cyclic loading 
to a constant value in the ratcheting regime. 



dissipated energy increases as Aa also increases. For very 
high loads, in fact, the dissipated energy can be even 
higher than the work exerted on the material (Aa — 0.45 
of Figure 7). This implies an overall modification of the 
structure, including of course the destruction and creation 
of contacts. Since all the curves have a plateau for a large 



enough number of cycles, a measure of the final value, Ap, 
can be worked out. This is shown in part (b) of Figure 
7. Above the shakedown limit, the increase of dissipated 
energy with Aa is potential in the range of parameters 
covered by our simulations. For very low Aa the system 
approaches the plastic shakedown regime, in which Ap is 
independent of Aa. If the applied load is such that no 
incremental collapse is observed in the response, the post- 
compaction regime is characterized by a growth of the 
accumulated strain with N as P(N) oc arctanipN + c) , 
being a, b, and c adjusting parameters. Figure 8 shows the 
fit of these data to the corresponding best-fit curves. This 
dependence is consequence of a simple second order poly- 
nomial relationship between the strain rate and the accu- 
mulation of strain. Only in case Aa — 0.40, some devia- 
tion is observed at the beginning of the loading, showing 
that this simple approximation is not valid for high loads. 
The behavior of the plastic strain in the post-compaction 
is related to a decay of the relative number of sliding con- 
tacts in the early phases of the loading. Due to the cyclic 
loading, a periodic sliding of the contacts is forced in the 
sample. The maximum number of sliding contacts in one 
cycle, N s , decays to a nearly constant value in the first 
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Fig. 8. Fit of the strain-rate vs. strain curves in the post- 
compaction phase to a simple law (see text), to some of the 
curves plotted in Figure 6. 



stages of the post-compaction. Figure 9 shows the behav- 
ior of N s scaled with the total number of contacts: n s . 
The final value of n s is clearly larger for the highest load- 
ings, and so is its variability. For Aa = 0.1 the number of 
sliding contacts remains practically unaltered during the 
ratcheting. For A = 0.4, however, this number changes 
continuously from one cycle to the next. This mainly re- 
flects the changes in the number of sliding contacts in the 
sample, although also the creation and destruction of con- 
tacts between particles is responsible for the changes in n s . 
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Fig. 9. Decay of the maximum relative number of sliding con- 
tacts n a in the post-compaction phase for the different loadings 
of Figure 5. 



The evolution of the resilient strain will also help to get 
an insight into the mechanisms involved in the deforma- 
tion. It is known in soils mechanics that the resilient strain 
changes only slightly through the simulation [10]. It can 
be inferred from Figure 10 that the adaptation time of the 
system depends on the imposed loading. For large load- 
ings, the resilient strain is still changing from one cycle to 



Yr (%) 



0.1 ■ 



A (7=0.32 
A a=0.36 

A o=0.4 

A o=0.5 * 



OOGOOGOOOOOOOOGOOOOOOOOOOOOOOOOOOOOOOGOOO 
A AA A AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 



50 



100 150 

N 



200 



250 



Fig. 10. Evolution of the resilient strain throughout the sim- 
ulation for some of the systems included in Figure 6. 



the next, after 200 cycles. For lower excitations, however, 
the resilient strain does not change after 20 or 30 cycles. 
In the ratcheting regime, the stress-strain cycle maintains 
its shape and the resilient modulus remains therefore con- 
stant. In the incremental collapse, however, the resilient 
strain varies, reflecting the changes in the configuration of 
the sample. The changes, being minor for moderate devi- 
atoric stresses (Aa — 0.4), become obvious when Aa is 
increased. 



Conclusions. 

A given material at rest or equilibrium under fixed stress 
conditions shows a certain configuration of contacts be- 
tween particles and forces. When a cyclic loading is im- 
posed, the system reacts by changing its configuration. In 
a continuous and gradual increase of the loading ampli- 
tude Aa, the material will start by trying to change the 
force skeleton without altering its configuration of con- 
tacts [19]. Negligible creation or destruction of contacts 
will occur and, for low enough excitations, no sliding con- 
tacts can be found in the simulation. After one cycle, 
however, the particles will recover its original state only 
approximately because some energy will have been dissi- 
pated, i.e. due to viscosity. At this range of loadings, the 
dissipated energy is independent of the loading and, basi- 
cally, does not change from one cycle to another. No pure 
elastic response is therefore found in the system, and the 
possibility of elastic shakedown is consequently also dis- 
carded. It is possible, however, that the material reaches a 
state in which there is no further systematic accumulation 
of permanent stress (plastic shakedown). 

For higher loadings, the energy input is first quickly 
dissipated by a re-arrangement of the sliding contacts of 
the system, the so-called post-compaction. The dissipated 
energy per cycle relaxes then to a stationary value depen- 
dent mainly on the loading, but also on characteristics of 
the grains such as the friction or the stiffness of the con- 
tacts. In this stage, the system evolves quasi-periodically. 
The sliding contacts show a periodic behavior within the 
cycles, although there is a linear increase of the deviatoric 
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strain and an overall compression of the system. This re- 
flects the fact that the system is still not dissipating all 
the energy it should. It is therefore expected that the ma- 
terial evolves on a much longer scale to a final shakedown 
state in which all the energy supplied to the system is dis- 
sipated. This process may take a longer time in the simu- 
lation than in the real experiment, where more dissipative 
mechanisms exist. 

If the loads imposed are high enough, there is no pos- 
sibility for the system to re-arrange to the new situation 
and the post-compaction is substituted by an incremen- 
tal collapse. The energy that the system is not able to 
dissipate in its configuration, is used to modify its shape. 

The stress-strain relations have been carefully studied 
in this paper, both as the mean to calculate the dissipated 
energy in each cycle, but also as the best way of charac- 
terizing the state of the system. In this respect, it has 
been proven that, for non-collapsing materials, the basic 
assumption expressed in (1), allows to fully characterize 
any stress state of the system in term of the strain rate, 
the resilient modulus and Poisson's ratio. Although there 
exist a large number of valuable models predicting the de- 
pendence of these resilient parameters on the stress [20], 
little is known about how they are affected by material 
constants such as friction or the stiffness of the contacts. 
This dependencies, however, are required in order to pre- 
dict the behavior of the UGM within a certain structure, 
and will be presented elsewhere. 
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